Static versus dynamic fluctuations in the one-dimensional extended Hubbard model 
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The extended Hubbard Hamiltonian is a widely accepted model for uncovering the effects of strong 
correlations on the phase diagram of low-dimensional systems, and a variety of theoretical techniques 
have been applied to it. In this paper the world- line quantum Monte Carlo method is used to study 
spin, charge, and bond order correlations of the one-dimensional extended Ifubbard model in the 
presence of coupling to the lattice. A static alternating lattice distortion (the ionic Ifubbard model) 
leads to enhanced charge density wave correlations at the expense of antiferromagnetic order. When 
the lattice degrees of freedom are dynamic (the Hubbard-Holstein model), we show that a similar 
effect occurs even though the charge asymmetry must arise spontaneously. Although the evolution of 
the total energy with lattice coupling is smooth, the individual components exhibit sharp crossovers 
at the phase boundaries. Finally, we observe a tendency for bond order in the region between the 
charge and spin density wave phases. 

PACS numbers; 71.10.Fd, 71.30.+h, 02.70.Uu 



I. INTRODUCTION 

The study of strong interaction effects in low- 
dimensional systems remains one of the most active fields 
of research in condensed matter physics. The extended 
Hubbard Hamiltonian (EHH) has been widely explored 
as a model of correlation effects in tight-binding systems 
and, more specifically, for the competition between dif- 
ferent types of ground state order: charge density wave, 
antiferromagnetism, and, in the case of attractive inter- 
actions, superconductivity. In one dimension, it has also 
been used to understand the behavior of materials includ- 
ing conducting polymers^ and organic superconductorsi^ 

The ground state phase diagram of the one- 
dimensional EHH was first obtained within a weak cou- 
pling renormalization group (RG) calculationi^i^ For re- 
pulsive on-site interactions U which are sufficiently large 
compared to the intersite repulsion V, specifically, for 
U > 2V, the ground state is a spin density wave (SDW) 
phase, with power law decay of spin correlations. For 
2V > U, the ground state has charge density wave 
(CDW) order. These charge correlations exhibit true 
long range order, that is, they go asymptotically to a 
nonzero value at large separations, since the associated 
broken symmetry is discrete. Finally, for attractive in- 
tersite interactions, singlet and triplet superconducting 
phases exist at T = 0, again with power law decays of 
the associated correlation functions. 

Subsequent to the RG work, the question of the or- 
der of the transitions between these different phases was 
studied, with a prediction that for repulsive U and V 
second-order SDW-CDW transitions at weak coupling 
were separated by a tricritical point from first-order tran- 
sitions at strong couphng.'^'^-''''^ Up to several years ago, 
estimates of the location of the tricritical point varied 
from Ut = IM to Ut = 5t (with Vt « Ut/2.) More 
recently, this picture has been further modified by the 
suggestion that a narrow region exhibiting "bond ordered 
wave" (BOW) correlations separates the SDW and CDW 



regions at weak coupling i^ i-'^^'^^i-'^^'-'^'^i-'^'^'^^ 

The competition of CDW and SDW order in the one- 
dimensional EHH is further modified if the electrons cou- 
ple to lattice degrees of freedom. In the case where these 
are static, most investigations have addressed the case 
when there is only on-site repulsion U, that is, ^ = 0. In 
this "ionic Hubbard model" the frozen distortions have 
an alternating pattern down the chain, and an addi- 
tional issue is the possibility that the band insulator at 
U = and half filling is first driven metallic before be- 
coming a SDW Mott insulator liiii^ If the coupling of 
the electrons to the lattice is in the form of dynamically 
varying phonon degrees of freedom, one has the Hubbard- 
Holstein or Su-Schrieffer-Heeger Hamiltonian. 

The interplay between band-insulating behavior and 
electron-electron interaction effects such as those stud- 
ied in this paper has recently been explored in a num- 
ber of contexts. Dynamical mean field theory studies of 
binary alloy band insulators described by a bimodal dis- 
tribution of randomly located one-body potentials have 
observed several novel effects, including Mott insulating 
behavior away from half fillin g^^i^° and band-insulator to 
metal transitions driven by increasing on-site repulsion^ 
Analogous studies of interacting bosons in "superlattice" 
potentials in which the site energies are modulated have 
also been used^Si'^'^iiS^i^^ to describe experiments on ul- 
tracold optically trapped (bosonic) atomsj ^^i^^'^^i'^° 

There has been relatively little work, especially using 
quantum Monte Carlo (QMC) simulations, which ad- 
dresses how such lattice coupling affects the SDW-CDW 
phase boundary in the EHH in which both U and V are 
nonzero. In this paper, we apply the world-line QMC 
(WLQMC) method to the one-dimensional EHH with an 
additional, static one-body potential, and with dynam- 
ically fluctuating ("Holstein") phonons. We quantita- 
tively determine the amount of lattice coupling required 
to stabilize a charge ordered phase when the system be- 
gins at values of the electron-electron interactions in the 
spin density wave regime. An interesting feature of our 
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results is that the quantum fluctuations induced by the 
hopping t have the opposite effect on the strong coupling 
{t = 0) phase boundary in the two cases. We also present 
detailed results for the evolution of the different compo- 
nents of the energy through the phase transition region. 

The remainder of this paper is organized as follows. 
An explicit description of our Hamiltonian and a brief 
review of our numerical approach are presented in Sec. 
im Results for coupling to static and dynamic lattice 
deformations are given in Sees. IIIII and llV| respectively. 



II. MODEL AND COMPUTATIONAL 
METHODS 

The extended Hubbard Hamiltonian is 
Hci = K + P, 

P ^ U'^Tii^^rii^l +V'^nini+i. (1) 

i i 

Here c|^, c^.o-, and ni^„ are the creation, destruction, and 
number operators, respectively, for electrons of spin a at 
site « of a one-dimensional lattice, and Ui — n^^^ -I- rii^i. 
The hopping t determines the kinetic energy (noninter- 
acting band dispersion ek — —2tcosk), and is set to 
t = 1. U and V, taken to be positive, are the on-site 
and intersite repulsions. We will be exclusively inter- 
ested in the properties of the model at half filling where 
the number of fermions Nf = = is equal to the 

number of lattice sites. 

We will consider additional couplings to an on-site lat- 
tice degree of freedom, 

H = Hcl + i?lattice, 

i/lHM = A^(-1)V„ 

i 

i i ^ ^ 

where i?iatticc can take one of two possible forms: either 
static (ionic Hubbard model "IHM") or dynamic ("Hol- 
stein" ) . Analytic and numeric studies on such Hamiltoni- 
ans are quite numerous ^^'^i^'^^'^s^^^^dodida 

It is useful to review the strong coupling {t — 0) phase 
diagram, since when the hopping is nonzero the topology 
of the phase diagram is rather similar qualitatively and 
even quantitatively. In the absence of an interaction with 
the lattice, the SDW phase, which consists of a collection 
of singly occupied sites, has energy i?s"5w ~ while 
the CDW phase has alternating empty and doubly occu- 
pied sites, and energy -Ecdw ~ The boundary is 
given hj V = U /2. A static lattice distortion A breaks 
the twofold symmetry of the CDW state and lowers the 



energy by A^A on the preferred sublattice. The resulting 
boundary is shifted to ^ = U /2 — A. 

In the case of coupling to a dynamical phonon, we can 
construct the t — Q phase diagram by completing the 
square of the electron-phonon term in the Hamiltonian. 
The result is an oscillator with the same frequency ujq and 
an equilibrium position shifted by A/wq. An attractive 
on-site interaction — (A^/2ajQ)ni jUi^i is also generated. 
Other terms can be absorbed into a shifted chemical po- 
tential and energy. As with the static term, the weak- 
ening of the on-site U shifts the strong coupling phase 
diagram in favor of CDW order. If — (A^/2a;Q) is suffi- 
ciently large, pairing correlations can come to dominate, 
especially in the doped case. We will not work in that 
parameter regime here. 

In order to understand how the quantum fluctuations, 
which develop as t increases, modify these simple con- 
siderations, we employ the world-line quantum Monte 
Carlo (WLQMC) method.'*'^ Consider first the approach 
for H = Hci + HiuM- We begin by discretizing the inverse 
temperature /? into intervals e — P/M in the partition 
function, and approximating the incremental (imaginary) 
time evolution operator by the product of the exponen- 
tials of the kinetic energy and potential energy terms 
separately. 

Z = Tr {e-"") « Tr (e-^^e-^(^+^™))^^. 

This Suzuki- Trotter approximation^^*^ introduces errors 
in measurements^^^ which are of order the commuta- 
tor [K,P], that is, tUe'^, tAe"^, and tVe'^. Except 
where otherwise noted, we will choose e = 0.25, which 
is sufficiently small that the systematic Trotter errors in 
the location of the phase boundary are comparable to 
those arising from statistical fiuctuations in the Monte 
Carlo sampling and uncertainties associated with finite 
size scaling. 

The construction of a path integral for Z is com- 
pleted by introducing complete sets of fermion occupa- 
tion number states / = l"-i.o-}r ("«,o'lr both for the 
trace and at all imaginary times, i.e., between each prod- 
uct, e^'^^e"^^^^^™"). The exponentials of the terms in 
P + HiHM immediately act on the eigenstates, replacing 
all operators by numbers. Thus the weight of a particu- 
lar occupation number configuration gets a contribution 
Wp Wmu, 

Wp{{ni^r,a}) = exp(^e^ [Uni^r,]'ni,T,i 

WiHM({?ii,r,dr}) = exp(^e^ A(-l)*(nj,7.j +"-i,r,i))) 

where {rii.r.o-} denotes the space- and imaginary-time- 
dependent occupation numbers in the collection of inter- 
mediate states. 



25 



20" ^ i 

() O ^ A A 



15 

[ 

10 
5 




N=8 □ 
N=16 o 
N=32 A 



0.5 



1.5 

M 



2.5 



200 
160 
120 

3 

80 
40 



1.2 1.3 
A/t 



1.8 
1.4 
1 

0.6' 



0.5 



1 



1.5 

A/t 



2.5 



22 



□ □ o 



10" ° 



6 « 



0.5 



1.5 

M 



2.5 



FIG. 1: Spin density wave (top), charge density wave (mid- 
dle), and bond ordered wave (bottom) susceptibilities versus 
staggered site energy A for = 6t, V = 1.5t, pt = 8, and 
= 8, 16, 32. The SDW-CDW transition occurs at close to 
the t = value, A = U /2 — V . BOW correlations are en- 
hanced in the intermediate region. In the inset to the central 
panel, the scaled xcdw is shown for 7 = 1. The scaled sus- 
ceptibilities cross at /S.c/t = 1.278, indicated by the vertical 
dotted line. 



To accomplish the same replacement of operators by 
numbers for the kinetic energy exponentials, K is further 
subdivided (the "checkerboard decomposition" ) 



I into 



K = Kodd + Kcven, 
i odd, cr 



I CVC11,CT 



The expectation value of Kodd and Kcvcn between the 
occupation number states \ni^a)^ and (ni^o-l^+i then re- 
duces to a product of independent two-site problems 
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FIG. 2: Illustration of a BOW phase. The thick lines indicate 
a high kinetic energy while the thin lines indicate a small 
kinetic energy. 



which can be solved analytically. Since particle num- 
ber is conserved in each hopping process, the number of 
electrons on each pair of sites in the two states to the 
left and to the right of the exponential is identical. Thus 
the world lines generated by connecting all occupied sites 
{ni^T,(y = 1) are continuous. The four nonzero matrix el- 
ements are 

(00|e'*("5"^+"2=i) |00) = 1, 
(11 |e^*('=I'=2+=2Ci) I 11) = 1, 
(10|e'*('=5'=^+"^=i) I 10) =cosh(te), 
( 10 I e^*('=l'=2+=2ci) I 01 ) = sinh(te). 

The product of all these factors over the space-time lat- 
tice constitutes a second contribution Wk to the weight 
associated with the configuration. Thus, the total weight 
is Wtot = WpWihmWk- Because all of the matrix ele- 
ments are positive in one dimension, the WLQMC algo- 
rithm does not exhibit a sign problem. 

In the case H — i?ci + -ffHoistoin, the trace and interme- 
diate states include not only fermion occupation labels, 
but also a complete set of phonon position eigenstates. 
As with iJoi, the exponential of the phonon kinetic and 
potential energies is discretized and split apart. The re- 
sult is that in addition to the electronic contributions 
Wp Wk there is a final phonon piece. 



W^ph({a;i,r}) = exp 



1 



2,r 



,2 2 



Let us then summarize the basic features of the sim- 
ulation. The degrees of freedom being summed over 
are two space-time arrays of occupation numbers rii^r,] 
and ni.T-.|, and, in the Holstein case, a space-time ar- 
ray of phonon coordinates xi^t, with i = l,2,...,N 
and T = 1,2,...,2M. (The factor of 2 comes from 
the checkerboard decomposition.) The total weight of 
the configuration is Wtot — WpWxWph- The elemental 
Monte Carlo moves consist of local distortions of the con- 
tinuous world lines, together with updates of the phonon 
degrees of freedom. Moves are accepted or rejected ac- 
cording to the Metropolis algorithm: a random number 
< r < 1 is generated and the move is accepted if 
r < WU/Wtot- 

The WLQMC algorithm can suffer from long autocor- 
relation times. Other approaches such as the stochas- 
tic series expansion metho d^^'^°i^^ and loop algorithms— 
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FIG. 3: Spin density wave (top), charge density wave (mid- 
dle), and bond ordered wave (bottom) susceptibilities versus 
staggered site energy A for U — 6t, V = O.Ot, {0.5t), 3.0t, and 
N = 32. 



can be used to speed up the evolution in phase space. 
Here we confine ourselves only to introducing global 
moves^ - in the phonon degrees of freedom to address even 
more serious large autocorrelation times there. 

We conclude with a discussion of the observables we 
will measure. The various components of the energy ex- 
hibit sharp features as the phase boundaries are crossed. 
Real space spin, charge (relative to the mean) , and bond 
operators are defined by 



n{l,T) = ni^r,i +nLT,i - 1, 

Hl,r)=y] (c|+, ^ (r)Q,, (r) -|- (r)Q+i,, (r)) . 



80 
60 

T ^ — 
40 



20 




o e o e G 



-40 

38 
34 
30 
26 
22 
11 

80 
60 
40 
20 



(a) 



>f....x-.-x.-.x-..x....x....^....><....^..-e::^ 



f Q H-B.-H-..H.^..H....a....g g 

20 t V V V ? V V ^ ^^ 








0.5 



1 



1.5 
A/t 



2.5 



^i....A 
St-*" 



,® n X X. I : B- 




"A 

ct> e e ie- 





v=o.ot - 


— ^ 


(c) 


V=0.5t 


■B- 




V=1.0t 


■■ X- 




V=1 .5t 


e 




V=2.0t 






V=2.5t 


...A-. 


.s... 


V=3.0t 




K "®\ 












■A... m. « 


K 





[3 □ H ■■B--B-B-E 



V ' V ' V ' ? ' V ' ^ ' V ' ? ' V ' V ' ^I? ' V ' V ' V ' V ' 







0.5 



1 



1.5 
A/t 



2.5 



(d) 




100 
80 
60 
40 
20 




FIG. 4: Total (a), kinetic (b), intersite Coulomb (c), and 
on-site Coulomb (d) energies at fixed U = 6t and different 
V. The inverse temperature is fixed at /3t = 8 and the lattice 
size at iV = 32. The kinetic energy is largest in the region 
where there is a balance between the CDW and SDW insu- 
lating tendencies, in good correspondence with the behavior 
of xbow. The intersite interaction energy falls abruptly on 
entry to the CDW state, and the on-site energies rise steeply 
as the pairs form. 
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FIG. 5: Phase diagram in the intersite V and staggered site 
energy A plane, with U = 6t and fUt = 8. Line with symbols 
is the result of the WLQMC simulations in this paper. We 
also show the exact result (line without symbols) for t — 0. 
As expected, the strong coupling limit works well at large V, 
but there are significant deviations as V becomes smaller. 



The associated correlation functions are 
Cspin(', t) = (m(/, T)m(0, 0)} , 

Cctiarge(',T) = {n{l , T)n{0, 0)) , 
^boiid 

where (0,0) is some reference site in our system. The 



local moment is defined as {ml) 



L-spin 



(0,0). 



We will also look at the Fourier transforms of these 
quantities. The equal time spin structure factor is 



*S'spin(9) 



1 



N ^ 

I 



spi 



„a,o). 



with analogous definitions for S'cimrgc and S'bond- 
corresponding zero-frequency susceptibility is 



The 



SDW susceptibility decreases and the CDW susceptibil- 
ity grows. Indeed, xcdw rises dramatically in the vicin- 
ity of A = [7/2 — as suggested by the strong couphng 
analysis. The transition becomes increasingly sharp as 
the lattice size is increased. Because the CDW correla- 
tions break a discrete symmetry, true long range order is 
possible at T = 0. With our normalization conventions 
we expect the CDW structure factor and susceptibility 
to grow linearly with lattice size after the onset of long 
range order. This is borne out in the central panel of 
Fig. [TJ The inset in this panel shows a scaled version of 
the raw data for xcdw- A crossing of the curves for dif- 
ferent lattice sizes N allows us to determine the location 
of the critical point. 
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again with analogous definitions for Xciiargc and Xbond- 

At half filling the largest responses in the Hubbard 
model are at wave vector q = n. In a disordered phase, 
c(l,0) decays exponentially to zero with the site separa- 
tion I, and the structure factor is independent of lattice 
size N. If true long range order develops, then the struc- 
ture factor grows linearly with lattice size, with the factor 
giTTi providing the necessary phases so that the oscillat- 
ing c{l) add constructively. The susceptibility similarly 
examines the asymptotics in imaginary time, diverging 
with P when c{l, r) remains nonzero for large r. 



III. RESULTS: EXTENDED IONIC HUBBARD 
HAMILTONIAN 
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FIG. 6: Spin density wave (top), charge density wave (mid- 
dle), and bond ordered wave (bottom) susceptibilities versus 
electron-phonon coupling A for U = 6t, V = 1.5t, uo = 2t, 
Pt = 8, and TV = 8, 16, 32. 



In the extended Hubbard Hamiltonian with U ^ 6t 
and V = 1.5t, we are well within the SDW phase since 
U > 2V. In Fig. [T] we see that as A is increased, the 



The SDW correlations that are dominant at small 
A break a continuous symmetry, and hence in one di- 
mension decay with a power law at T = 0, that is. 
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Cspin(',0) cx l/l. This behavior accounts for the rela- 
tively less rapid growth of the SDW susceptibility with 
lattice size. 

It is important to make another distinction between 
the CDW and SDW phases, the phases that arise as bro- 
ken symmetries from the interaction terms V and [/, and 
the staggered density which is caused by the one-body 
term A. This staggered potential A breaks translational 
invariance so that there is a small degree of CDW order 
even in the SDW phase. By contrast, in a competition 
solely between U and 1/ at A = 0, no CDW order would 
exist in the SDW phase. 

The bottom panel of Fig. [T] shows the BOW correla- 
tions. In a BOW phase the kinetic energy on the links 
oscillates between two values as one traverses the chain 
(see Fig. [2]). SDW correlations are immediately plausi- 
ble after observing that U leads to singly occupied sites 
(moment formation) and that neighboring spins that are 
antiparallel have a second-order lowering of their energy 
{/S.E'^'^^ cx —t^/U) relative to neighboring spins that are 
parallel. Analogous reasoning applies to CDW correla- 
tions. A picture of the less familiar BOW order is the 
following: consider a CDW pattern of doubly occupied 
and empty sites. A fermion hopping from doubly occu- 
pied site i onto neighboring empty site will prevent, 
through the interaction U, the hopping of a second elec- 
tron from doubly occupied site i -\- 2 onto Instead, 
an electron on site i + 2 would prefer to hop to i-|-3. Thus 
the bonds (i,i -I- 1) and (i -|- 2, i -I- 3) have high kinetic 
energy, while the intermediate bond {i+l,i-\-2) has low 
kinetic energy. This way of understanding the origin of 
BO invokes both CDW and SDW correlations, making 
it plausible that the BOW might form on the boundary 
between the two. 

Since the BOW phase also breaks a discrete trans- 
lational symmetry, the associated ground state order 
should be long ranged. As mentioned in the Introduc- 
tion, in the extended Hubbard model (A = 0) the orig- 
inal picture of the phase diagram was one with only 
SDW and CDW regions, with a weak coupling second- 
order transition changing at a tricritical point to a strong 
coupling first-order transitioni^i^i^ Recent QMC simu- 
lations with the stochastic series expansion (SSE) have 
suggested instead that, at weak coupling, as V is in- 
creased at fixed U there are two separate transitions: 
a SDW-BOW transition of the Kosterlitz-Thouless type, 
followed by a second-order BOW-CDW transition. These 
transitions merge at a multicritical point into a single, 
direct, first order SDW-CDW transition line at strong 
coupling. The multicritical point was found to be at 
{Um,Vm) = (4.7 ±0.1, 2.51 ±0.04). 

Aspects of this conclusion had been challenged by den- 
sity matrix renormalization group calculation a^^i^^'^^ . In 
particular, the suggestion is that the BOW phase exists 
only precisely on the SDW-CDW transition line, as op- 
posed to being present in an extended region. Moreover, 
rather than starting a± U = V = Q and reaching out to 
the multi-critical point, the BOW line was concluded to 
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FIG. 7: Spin density wave (top), charge density wave (mid- 
dle), and bond ordered wave (bottom) susceptibilities versus 
electron-phonon coupling A for U = 6t, V = O.Ot, {0.5t),3.0t, 
LOO = 2t, and = 32. 



begin at finite, nonzero coupling and also extend some- 
what beyond the numerical value for the multi-critical 
point obtained using the SSE. Further SSE calculations ^2, 
and functional RG treatment&l^ appear to confirm earlier 
SSE work. 

We do not propose here to add to this discussion, since 
our main focus is on the shift in the SDW-CDW phase 
boundary. Indeed, the value of U in Fig.[T]is large enough 
that we would likely be above the BOW region of the 
phase diagram. Nevertheless, the bottom panel of Fig. [T] 
does indicate a pronounced maximum in xbow near the 
SDW-CDW transition, hinting that such order may be 
present at weaker coupling. If long range BO were to 
exist, we would expect to see xbow grow linearly with 
N, as does xcdw- This is clearly not the case for the 
parameters and lattice sizes of Fig. [TJ 

In Fig. [3] we fix J7 = 6t and the lattice size at — 32, 
and sweep A for different choices of V. As expected, the 
size of A required to destroy the SDW phase decreases 



7 




-40 ~ 

1 2 3 4 5 6 

U\ 

FIG. 8: (a) Kinetic (hopping), (b) intersite Coulomb, (c) on-site Coulomb, (d) kinetic (phonon), (e) potential (phonon), (f) 
electron-phonon coupling, and (g) total energies versus coupling constant A for U — 6t, V = O.Ot, (0.5t),3.0t, pt — 8, ujq = 2t, 
and A'' = 32. The fluctuations in the phonon kinetic energy are significantly smaller than those for all other energies for these 
parameters. 



as the intersite interaction V , which cooperates with A, 
rises. As with the data of Fig. [U the fall of xsdw co- 
incides closely with the rise of xcdw- In each case the 
transition is marked also by a maximum in xbow- The 
sharpness of the peak in xbow diminishes as V grows, 
which is consistent with the SSE'^^ and density matrix 
(DM) RG^ calculations on the extended Hubbard model 
which (although they disagree in certain respects) both 
conclude that BO is not present at strong coupling. We 



note that a Mott-insulator-BO transition has also been 
suggested by Zhang et al. in the V = {) limit with A = 2.0 
and Ua = 5.95±0.0liS 

The behavior of the total energy, Fig.dKa), is feature- 
less through the sweep upward in A. However, abrupt 
evolution of the individual components of the energy, 
Figs.[3]Jb)-(d), accompanies the transitions in the suscep- 
tibilities. The energy associated with V decreases sharply 
upon exiting the SDW phase where adjacent sites are oc- 
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FIG. 9: This figure shows the phase diagram for intersite V 
and electron-phonon coupling A with U — 6t, cuo ~ and 
/3t = 8. Line with symbols is the result of the WLQMC 
simulations in this paper. The functional form of the fit is 
V = + b. We show the exact result (line without symbols) 
for the t — phase. 



cupied, while the energy associated with U jumps upward 
with the development of double occupancy. The kinetic 
energy is relatively benign, but, like xboWi reaches max- 
ima along the SDW-CDW transition line. Evidently, the 
near balance between the insulating tendencies of U and 
V allows greater fluctuation in the electron positions. 

The values of A at which the different susceptibilities 
change abruptly, and at which features in the energy are 
also evident, enable us to draw the phase diagram in the 
V-A plane for fixed U = 6t shown in Fig. O At A = 
our QMC results match quite nicely the DMRG results 
of Jeckelmann^ who finds 14 = 3.155±0.005 for U = 6t. 
This A = transition point is not too far shifted from 
the strong coupling value Vc = U/2 = 3t when U — 6t. 

As the staggered potential A becomes greater, our 
QMC phase boundary bends more away from the t — 
line Vc — U/2 — A. The SDW phase appears to terminate 
at A = 2.46±0.05 in the absence of intersite repulsion V. 
While labeled as CDW, the large A phase in the V = 
limit is perhaps more properly termed a band insulator, 
where the alternating charge density is a consequence of 
the staggered one-body potential as opposed to many- 
body effects. 



IV. RESULTS: EXTENDED HUBBARD 
HOLSTEIN HAMILTONIAN 

Having completed our discussion of the case of the in- 
terplay of a static alternating external potential with the 
correlation terms U, V in the extended Hubbard Hamilto- 
nian, we now give analogous results for the case when we 
couple to dynamical (Holstein) lattice distortions. Fig- 
ure [6] is a companion to Fig.[Tl showing the evolution of 
the spin, charge, and bond susceptibilities with electron- 
phonon coupling A (rather than staggered potential A) 
for different system sizes N. As discussed earlier, A has 



a similar qualitative effect to A, since it weakens the 
on-site repulsion U and hence drives CDW formation. 
There are significant quantitative differences between the 
two situations. The SDW-CDW transition as a function 
of electron-phonon coupling A appears to be much more 
abrupt. Recall that A breaks the lattice symmetry ex- 
plicitly, selecting out a single preferred sublattice. It in- 
duces CDW order even within the SDW phase and as 
a consequence the change through the transition is less 
dramatic. The Holstein interaction, in contrast, sponta- 
neously breaks the translational symmetry when it drives 
CDW order. We note further that BOW order is less 
sharply peaked at the SDW-CDW boundary. 

Figure[7]is a companion to Fig.[3l similarly showing the 
susceptibilities as a function of electron-phonon coupling 
constant A for a collection of values of ^ at a single lat- 
tice size TV = 32 and ljq = 2t. As V increases, a smaller 
A is sufficient to drive CDW formation. There appears 
to be some variation of the sharpness of the evolution of 
the susceptibilities near Ac as V is varied, with the most 
abrupt behavior occurring for intermediate V . In the ex- 
tended Hubbard model (A = 0), the transitions become 
monotonically more steep with increasing V. Indeed, as 
noted earlier, they change from continuous to discontin- 
uous beyond the tri- (multi)critical point. The fluctua- 
tions of XCDW at large A in Fig. [7] (middle panel) often 
occur in QMC studies of electron-phonon Hamiltonians 
and are associated with long equilibration times which 
occur when the electrons and lattice degrees of freedom 
are strongly coupled. 

As with HinM, the components of the energy (Fig. [8]) 
lend important supporting evidence for the locations of 
the transition points. The behavior of Ey and Eu is the 
same as that observed previously in Fig.[4l and is more or 
less clear: in the SDW phase most sites are singly occu- 
pied and there is a significant contribution to Ey, which 
then drops abruptly in the CDW phase where doubly 
occupied and empty sites alternate. In contrast, Eu is 
small in the SDW phase since sites are singly occupied, 
but then increases sharply in the CDW phase. What is 
perhaps less intuitive is the evolution of the phonon con- 
tributions to the energy. As A grows, the t — analysis 
suggests a smooth quadratic increase, £^phonon = A^/2cjg. 
Instead the phonon potential energy remains relatively 
flat throughout the SDW region, and then jumps up as 
the CDW is entered. The phonon kinetic energy is espe- 
cially interesting, showing a well-defined minimum in the 
transition region. The origin of this effect is not clear, 
-^phonon measured by the fluctuations of the phonon 
coordinates in imaginary time. Naively, one might ex- 
pect kinetic lattice fluctuations to be largest in the SDW- 
CDW transition region where the system is undecided 
between which type of order to assume. In the case of 
the electron kinetic energy we see precisely this effect in 
Fig. [HKa). The opposite appears to be the case for the 
phonon kinetic energy. 

Finally, Fig. [9] shows the phase diagram in the V-X 
plane at fixed U = 6t. It shares the same general fea- 
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tures as Fig. [5] with a SDW phase near the origin that 
is destroyed when either the intersite repulsion V or the 
electron phonon coupling A increases sufficiently. Figure 
ini describes how large a value of electron-phonon coupling 
A is required to convert the SDW phase, favored by U, 
to the CDW phase, favored by V, and is representative 
of how A affects the extended Hubbard model phase di- 
agram at all intermediate to large interaction strengths. 
It is important to note that, unlike Fig. [5l the QMC 
phase boundary does not bend away from the t — line 
Vc — U/2 — /4ojo- Instead, the boundary is uniformly 
shifted to increase the critical intersite repulsion, favoring 
SDW order. Again, the A = point on our phase bound- 
ary {Vc = 3.124 ± 0.011) agrees well with Jeckelmann's 
DMRG treatment. (See above discussion of Fig.[5l) 

For finite A we can compare against the phase dia- 
gram of Sil and Bhattacharyya who study the same ex- 
tended Hubbard model coupled to Holstein phonons4^ 
They draw the phase boundary in the U-V plane for 
different electron-phonon couplings. Translating to the 
units used in our paper, for U = 6t and V = 2t, their data 
suggest that the CDW phase is destroyed at Ac ~ 2.8. 
Our Fig. [HI gives Ac « 3.0 at V = 2t for the same pa- 
rameters. Likewise, Sil and Bhattacharyya find that for 
A = 5.6 there is no SDW phase at U — 3t. This is again 
nicely consistent with our data, which suggest that when 
A = 5.04 ± 0.06 there is CDW order. 



V. SUMMARY 

In this paper we have presented world-line quantum 
Monte Carlo simulations of the extended one-dimensional 



Hubbard Hamiltonian to which coupling to static stag- 
gered (ionic Hubbard) or dynamic (Holstein) lattice de- 
grees of freedom is added. The evolution of the suscepti- 
bilities to different types of order and the components of 
the energy were examined. For both static and dynamic 
couplings the region of charge density wave order in the 
phase diagram is found to be stabilized, and the phase 
boundaries are pinned down. Bond order is shown to 
be enhanced in the vicinity of the spin density to charge 
density transition. 

The results obtained in this work also show good 
agreement with previous studies. The zero-coupling 
limit (A, A = 0) results conform well with Jeckelmann's 
DMRG results. For dynamic couplings, the results com- 
pare favorably with the results of Sil et al. 

A comparison of the QMC phase boundaries with their 
counterpart in the t = limit shows a markedly dif- 
ferent behavior between the two types of coupling. For 
static couplings, the CDW phase is enhanced in the QMC 
calculation. Conversely, there is an enhancement of the 
SDW correlations with Holstein phonons. 
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